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A quantum phase transition is typically induced by tuning an external parameter that appears 
as a coupling constant in the Hamiltonian. Another route is to vary the global symmetry of the 
system, generalizing, e.g., SU(2) to SU(A^). In that case, however, the discrete nature of the con- 
trol parameter prevents one from identifying and characterizing the transition. We show how this 
limitation can be overcome for the SU(A'^) Heisenberg model with the help of a singlet projector 
algorithm that can treat A'' continuously. On the square lattice, we find a direct, continuous phase 
transition between Neel-ordered and crystalline bond-ordered phases at Nc — 4.57(5) with critical 
exponents z = 1 and jS/v = 0.81(3). 



I. INTRODUCTION 

The field of quantum magnetism encompasses a large 
variety of physical phenomena that are of current exper- 
imental and theoretical interest. These include competi- 
tion between interactions (frustration), ordering in con- 
ventional or unconventional magnetic states, and the ex- 
istence of fractionalized excitations. In two dimensions, 
where some of the most unusual physics occurs, there 
is a conspicuous absence of methods for studying the be- 
haviour of quantum magnets with high precision. On the 
analytical side, neither the powerful methods devised for 
dimension d = \ (bosonization, conformal field theory) 
nor the mean-field methods exact in high d are available. 
On the numerical side, simulations are difficult because of 
the enormous size of the Hilbert space and, for stochastic 
methods, because of the fatal "sign problem." 

One way to relax these strong methodological con- 
straints is to decrease the role of quantum ffuctuations. 
For instance, considering the classical limit of magnets 
with large spin S eases analytical studies. This limit, 
however, very often misses the important competition 
between the instabilities existing only at the quantum 
level. An alternative route — one that preserves the quan- 
tum fluctuations — is to enlarge the symmetry of the 
model, e.g., by extending the SU(2) spin symmetry to 
SU(A^). This has proved very useful in the past, as the 
N ^ oo limit often allows for an exact analytical treat- 
ment. Methods to study corrections are also avail- 
able, although they cannot fully capture the exact details 
of what happens at finite N. The advantage of SU(Af) 
models over classical ones is that they naturally allow for 
quantum states of matter (such as valence bond solids) 
by construction, even if the off-diagonal elements of the 
Hamiltonian are suppressed in the large- A?^ limit. 

Using this technique, and building on previous work,'^' 
Read and Sachde\0have studied in detail the SU(A^) gen- 
eralization of the Heisenberg Hamiltonian on the square 



lattice. For sufficiently large A^, the system sponta- 
neously breaks lattice translation symmetry to form a 
valence bond crystal (VBC). For small A^ [including the 
standard SU(2) model at A^ = 2], the ground-state is 
antiferromagnetically ordered. The details of the phase 
diagram and of the VBC depend on the representation 
of the generators of the SU(A^) algebra considered: for 
the case of square Young tableaux with n columns, a di- 
rect phase transition between the Neel and VBC states is 
predicted to occur at the (mean-field) value N/n ^ 5.26. 

In a technical breakthrough, Kawashima and cowork- 
e]- j3|4|5 | extended a quantum Monte Carlo (QMC) loop 
algorithm designed for SU(2) models to the SU(A^) case 
(for all integer A^ and for all single-row Young tableaux). 
Studying the square lattice case with this exact numeri- 
cal method, they found for n = 1 that the A^ = 4 model 
is Neel ordered, whereas the N = 5 model supports VBC 
order. This confirmed the analytical large- A^ predictions, 
but because of the discrete nature of the algorithm, these 
studies could not rule out an intermediate phase between 
N = A and 5. Even if this (possibly spin-liquicF') phase 
does not exist, it is impossible to obtain a precise value 
for the critical parameter A"c separating the two phases 
and to ascertain the nature of the phase transition at A'c. 

In this paper, we describe a quantum Monte Carlo al- 
gorithm, formulated in the total singlet basis, that can 
treat the parameter A'^ continuously (in the manner of 
analytical, large- A^ techniques). Applying this approach 
to the square-lattice SU(Af) Heisenberg model, we find 
that there is a direct transition occurring at A"c = 4.57(5) 
between the Neel and VBC columnar phases. The tran- 
sition is found to be second order, with critical exponents 
z ~ 1 and P/v ~ 0.81(3). At the end of the paper we dis- 
cuss the implications of finding a second-order quantum 
phase transition between states with incompatible sym- 
metries and the possible connection to the deconfined 
quantum criticality (DQC) scenario.'^ 



II. MODEL HAMILTONIAN 

Our starting point is the SU(-/V) generalization of the 
quantum Heisenberg model: 

J ^ 

h = -jJ2h,,^-J2J2 m^J^ij)- (1) 

a, [3=1 

Here, the exchange coupling J — 1 sets the energy scale, 
and denotes nearest-neighbor sites i and j. The 

generators of the SU(A^) algebra, J'^, satisfy the anti- 
commutation relation 

[Jj^i^), 4U)] = - <5„-^ J^",(z)). (2) 

We consider the "quark-antiquark" model, taking the 
fundamental representation of the generator on one sub- 
lattice (a single-box Young tableau) and its conjugate 
{N — 1 boxes in one column) on the other. The fusion 
rule for the two representations is depicted graphically in 
Fig-Sa). 

If we describe the fundamental representation using a 
basis |a) with a = 1,2, . . . ,iV, then the conjugate rep- 
resentation has states \a) that are fully antisymmetrized 
tensor products of the form 



\a2, 



. , aN) 



(3) 



An SU(A^) singlet formed between spins at sites i and j 
in opposite sublattices is the maximally entangled state 



1 ^ 



a I a 



(4) 
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Equation ([T|) can be understood as an operator that per- 
forms a local singlet projection Hij = —jfJ^{i)Ji{j) = 
\(f))ij{(t>\ij across all links of the square lattice. 

The SU(A^) Heisenberg model can alternatively be 
written as an SU(2) system with spin S = (iV — 1) /2 mo- 
ments interacting via higher-order exchange processes. 
An exact mapping connects the conventional SU(2) spin 
operators to the SU(A^) generators as follows: 
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The Hamiltonian can then be expressed in terms of 
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FIG. 1: (Color online) (a) An A-sublattice and a B-sublattice 
spin can pair to form a singlet, which corresponds to a column 
of zero (modulo A'^) boxes, (b) Update rules for the action of 
Hij (red line) on VB singlet states (black bonds) . 



For instance. 



Hi j - 



1/4 - S, • for S = 1/2 {N = 2); 
H,j = i[(Si • Sjp - 1] for = 1 (iV = 3); etc. This was 
the starting point of previ ous finite-temperature QMC in- 
vestigations of this model,!^ where a path-integral tech- 
nique was developed in the Sz basis of the spins. ^ In this 
paper, we take a rather different route, using a T = 
algorithm formulated in the SU(2) total singlet basis of 
the spins S. 



III. BOND BASIS 

Consider the subspace formed by bipartite valence bond 
(VB) stated in which two spins S in opposite sublattices 
form a VB by coupling pairwise in a singlet, formally 
given in the Sz basis by 



^/2S+l 



E (-1)' 



(9) 



[cf. Eq. ([4|]. For general S, this subspace does not span 
the full SU(2) singlet manifold, but only the subspace of 
states that are also SU(2S' -I- 1) symmetric. 

For bipartite valence bonds, we can impose a VB ori- 
entation convention such that the overlap between any 
two states is positive. This basis is nonorthogonal, and 
the overlap between two VB states is 



{vi\v2) - {2S+1) 



(10) 



where Ni is the number of loops formed by superimposing 
the two VB states \vi) and \v2), and is the number 
of VBs in each state (or, equivalently, half the number 
of lattice sites). This is a simple generalization of the 
well-known overlap rule for S = 1/2. 
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FIG. 2: (Color online) Loop diagrams contributing to two- 
and four-spin correlators [Si • Sj and (Si • Sj)(Sfc • S;), respec- 
tively] and their contributions. These overlap loops (schemat- 
ically circular here) are obtained by superimposing the VB 
configurations \vi) and \v2)- The straight (red) lines link the 
corresponding sites of the correlation function. Only diagrams 
with nonvanishing contributions are displayed. 



The Perron- Frobenius theorem tells us that on a bipar- 
tite, finite lattice the SU(iV) symmetric Hamiltonian (fl] 
admits a unique ground state. This state is an S\]{N ^ 
singlet, which can be expressed in the bipartite VB basis 



The operator Hij obeys the rules Hi 



and 



ijw/iiwikj = j4\(l>)iM)ki- As a consequence, the ac- 
tion of Hij on VB states is extremely simpl^and consists 
of the bond rearrangements depicted in Fig.[TJb). We ex- 
ploit this fact to simulate the SU(A^) model, noting that 
the VB projector QMC (Ref. ^ developed for 5 = 1/2 
works with precisely such update rules. For the sake of 
completeness, we describe this method (emphasizing the 



few details that differ) in Sec. IV 



What remains is to determine how to compute observ- 
ables. It is well-known in the S = 1/2 case that most 
observables can be written in terms of estimators based 
on the overlap loops.^^ For instance, the spin correlator 
(fijSi • Sj\v2) / {vi\v2) = (3/4)eij if spins i and j belong 
to the same loop, otherwise. Here e^^ = 1 if i and j 
are on the same sublattice, —1 otherwise. We have gen- 
eralized these rules for the spin S case, and the resulting 
nonvanishing loop diagrams are given in Fig.|2] alongside 
the value of their contribution. 

We now make our key observation: since the update 
rules, the Monte Carlo weights, and the estimators have 
analytical expressions in N (or S, via N = 25-1-1), 
simulations can be performed for arbitrary, continuous 
values of N. Formally, this can be understood as an 
analytic continuation from integer to real N. This is a 
great advantage over other QMC techniques,'^ which are 
restricted to half-integer and integer S. Large- analyt- 



ical techniqued^ can also treat continuous values of A'', 
but our numerical technique allows for an exact treat- 
ment of the Hamiltonian, in contrast to the mean-field 
approximation inherent to the large- approach. 

IV. NUMERICAL METHOD 

A. Quantum Monte Carlo algorithm 

We employ an efficient Monte Carlo algorithm that 
is an extension of the VB projector scheme introduced 
by Sandvik.flQ' The idea is to sample the ground state 
via the power method by applying (— _ff)*^ (for fixed M 
sufficiently large) to an arbitrary valence bond trial state: 



lim (- 

Af^oo 



(11) 



For the purpose of sketching out the algorithm, it is con- 
venient to rewrite Eq. ([T|) in a form that explicitly indexes 
the bonds that are acted upon: 



^h/j = J2h^, = J2h, 



(12) 



We make the equivalence Hh — Hi(^i,j j(^i,^ 



where b labels 

all the nearest-neighbour bonds on the square lattice. We 
now expand the powers of the Hamiltonian as 



i-Hr=(j2H, 



M 
b] = 



M 



{6i,b2,...,bA/}fc=l 



(13) 



Each sequence {61, 62, ^a/ } corresponds to the process 
in which \ipt) is acted on by local singlet projectors to 
give the propagated VB state \(pAi): 



H, 



.Hb^Hb^\ipT) 



(14) 



The sum over all possible sequences in Eq. ( 13 1 is eval- 
uated stochastically. For the SU(7V) model, the weight 
of each sequence is simply the product of the 1/A^ fac- 
tors (with iV = 2 for 5 = 1 /2) that appear as a result of 
the bond rearrangements [see Fig. [l]^b)] induced by the 
Hb's. The most basic Monte Carlo move then consists 
of replacing a few H},'s at random. Such changes are ac- 
cepted or rejected depending on the ratio of the weight 
before and after the move.'^ 

To compute observables, two such sequences are gener- 
ated with an additional factor in the sampling weight cor- 
responding to the overlap of the two propagated states. 
For the SU(7V) model, this is computed according to 
Eq. ( 10 1 The observables can then be measured using 
the loop estimators of Fig. |2] (see Sec. IV B I. 

The spin gap Ag, the energy difference between the 
first triplet excited state and the singlet ground state, 
can be computed using the same triple t pro pagation tech- 
nique introduced for the 5 = 1/2 case.l^^' Working with 
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a single propagation sequence, we reinterpret the initial 
VB trial state as containing one triplet bond; the only 
difference in the propagation rules is that the triplet is 
annihilated by Hh when acted on directly [i.e., the co- 
efficient 1 in the first rule of Fig. [ijb) is replaced by 0]. 
Looking at the statistics of the states that are not annihi- 
lated during the full propagation, one can easily compute 
the spin gap Ag. We refer the reader to Refs. 8 and ITUl 
for more details. 

To accelerate the convergence of the algorithm with 
respect to M, it is useful to sample a valence bond trial 
state that is a superposition of all valence bond configu- 
rations V with amplitude ipT^v): 



(15) 



There is no restriction on the trial state other than 
that (fT{v) be real and nonnegative and that ratios 
'^t{v2) / 'frivi) be easy to compute for small changes in 
configuration vi V2- In this work, we make use of a 
simple RVB trial stat^^^' in which the weight (Pt{v) — 



n 



hij is a product of individual bond amplitudes 



hij = 1 / rfj that fall off algebraically as a function of the 
bond length. The trial state is thus characterized by a 
single exponent p, which is a free parameter in our sim- 
ulations. To ensure that our results are fully converged 
(i.e., that M has been chosen sufficiently large) and do 
not show any residual dependence on the choice of p, 
we have carried out the ground state projection starting 
from four different RVB trial states (two magnetically 
ordered and two disordered) corresponding to p — 2.7, 
3.0, 3.5, and 5.0, and have checked that all observables 
converge to the same values. 

We should point out that the updates are simple and 
sign-problem- free only because Hij, appearing in Eqs. ([T]) 
and (|8]), is a pure, local-singlet projector. Our algorithm 
does not apply to general spin-S* Hamiltonians. In any 
case, the bipartite VB states do not form a basis for the 
singlet sector of general spin-S" models. 



B. Calculation of observables 

The quantities of interest at the transition are the stag- 
gered magnetization 



M=-l5](-l)'-^+'-='Sr 



its Binder cumulant 



U =1- 



3(M4) 
5(M2)2' 



(16) 



(17) 



and the dimer order parameter D = (D^, Dy), whose two 
components 

= Z^E(-l)''°Sr-Sr+e. (a = X , y) (18) 



are directed along the square lattice vectors and e^^. 

The measurements (M^), U, and (D^) are carried out 
using the two- and four-point rules shown in Fig. [2] For 
instance, using the same techniques as in Ref. ^ we ob- 
tain the following estimators: 



{^M\M^\ip'M) 

{'^mWm) 

and 



S{S+\)Y,Ll 



(19) 



{"PmWm) 



-[S{S+\)+2S\S+\f]Y,Ll 



-S{S+\)Y^Ll 

OL 



(20) 



where the sum on a runs over all loops formed by the 
overlap of the two propagated states, and Lq, 

is the size of a given loop. 



V. RESULTS 

We proceed by presenting our results for the square 
lattice SU( A^) model. Figure [s] shows the square of the 
staggered magnetization and its Binder cumulant. It is 
apparent that in the thermodynamic limit a dome of an- 
tiferromagnetic order survives for 1 < iV < iVc with 
between 4 and 5. The Binder cumulant U vanishes on 
the large-A^ side of the transition and appears to have 
a crossing point at N k, 4.3, which drifts rightward as 
L increases. We do not have data on sufficiently large 
systems to ascertain that the crossing is stable and can 
be unambiguously identified as the critical point. As we 
will see, a finite-size scaling analysis provides a better 
estimate of N^. 

The destruction of the Neel order is driven by the grad- 
ual elimination of singlet pairs that are correlated over 
long distances. In the large-A^ phase, the singlets are 
predominantly short-ranged and form domains of colum- 
nar ordering (see Fig. |4|. On small lattices, these or- 
dered domains are weak and they are almost equally dis- 
tributed among the four degenerate configurations. This 
can be seen in the ring-like structure of the probability 
distribution P(D,j,, Dy), shown in Fig.l5| that appears for 
N > Nc- Previous work^ * established the existence of 
VBC order for A'^ > 5, but could not determine whether 
it was of columnar [(D) ~ (0, ±D), (±1?, 0)] or plaque- 
tte [(D) ^{±D,±D)] symmetry. Our results suggest the 
former. There does not appear to be a true U(l) degen- 
eracy, as suggested in Ref. |3] Instead, it seems that there 
is an additional length scale ^vbc much larger than the 
spin correlation length ^. For ^ < L < ^vbC; there is 
merely an effective U(l) degeneracy. 
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FIG. 3: (Top-left panel) Square of the staggered magnetiza- 
tion M = 1)'^°' Sr as a function of for systems 
up to linear size L — 48 (lines are guides to the eyes). (Top- 
right panel) The Binder cumulant U measures the kurtosis of 
the staggered magnetization with respect to a purely Gaussian 
distribution. It vanishes in the limit L — > cx3 in the absence of 
antiferromagnetic order. (Bottom panel) A histogram of the 
magnetic observable is shown for values of spanning the 
transition. The distribution has a single peak. 
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FIG. 4: (Color online) Snapshots of typical VB configura- 
tions for a system of size L = 32. Short, nearest-neighbour 
bonds are drawn with a line; long bonds are indicated by cir- 
cles (whose area is proportional to bond length) at their end 
points. For N < Nc, many long bonds stretch across the sys- 
tem. For N > Nc, short bonds dominate. The shaded (pink) 
area marks a crystal domain with columnar bond order. 




iV=5.5 A^=6 




FIG. 5: (Color online) Density plots of histograms P{Dx, Dy) 
of the X- and y- components of the dimer order parameter 
(L = 32). On the magnetic side of the transition (e.g., A'' = 
4.5 < A^c), the distribution is characterized by a central peak. 
On the VBC side, the distribution is ring-like but develops 
additional weight along the main axes as N is increased. 



A finite-size scaling analysis of the magnetization and 
dimer order parameters suggests a continuous transition 
defined by a single critical value and a single set of 
critical exponents (see Fig. |6|. The data are not suffi- 
ciently sensitive to fix the exponent precisely, and the 
unusual behaviour of the Binder cumulant — its negative 
region and strong subleading corrections — makes it un- 
reliable for obtaining an independent estimate of v. Rea- 
sonable fits seem to be achievable for a range of values 
0.75 < v < 1. On the other hand, the ratio P/v is rel- 
atively stable. Repeating our fitting procedure for 
and independently with 3000 bootstrap samples, we 
conclude that both quantities vanish simultaneously and 
continuously at = 4.57(5) with P/iy = 0.81(3). Note 
that the anomalous dimension ?/ = 2/9/i/ — 1 = 0.63(4) 
is at least an order of magnitude larger than what would 
be expected from either the three-dimensional 0(3) or 
Z4 universality classes. 

In Fig. [t] we plot the spin gap as a function of 1/i 
in the region 4 < < 5. For A^ large enough, the data 
converges to a nonzero value in the thermodynamic limit. 
On the other hand, for the smallest A^, the gap clearly 
vanishes. A linear extrapolation in 1/L (dashed lines in 
Fig. [7]) reveals that the gap closes at A^ « 4.6, in good 
agreement with the value Nc = 4.57(5) derived from the 
order parameters. An unconstrained power-law fit yields 
a slope 1.02(10), confirming that the dynamical critical 
exponent of the phase transition is z = 1 . 

As is the case for any finite-lattice simulation, we can- 
not completely rule out an extremely weak first-order 
transition. Nonetheless, we have taken great care to ex- 



6 




-4 -2 2 4 4.4 4.5 4.6 4.7 4.8 

FIG. 6: (Color online) Simultaneous data collapse of the 
Neel and dimer correlations can be achieved with a single 
set of critical exponents. The left panel shows the QMC 
measurements plotted in rescaled coordinates with the val- 
ues V = 0.88, /3 = 0.71, and A''c = 4.57(5). Other values 
oi P/u ~ 0.8 produce good data collapse. The right panel 
shows bootstrapped results oi P/v versus A'^c when fits of 
and are performed independently. A single critical point 
corresponds to the crossing P/v = 0.81(3) and A^c = 4.57(5). 



tion or of a distribution with a complicated multi-peaked 
structure. The histograms discussed above rule out the 
latter: the negative U seems merely to correspond to a 
region in which the distribution of the magnetic order 
parameter is super-Gaussian (excess kurtosis). Another 
important observation relates to how the negative region 
evolves with system size. For first-order phase transi- 
tions, one typically observes very large negative values 
of the Binder cumulant that increase in magnitude with 
the system size. This has been explained phenomeno- 
logically by Vollmayr et aL^^^l For the transition observed 
here, the opposite is true: not only are the values only 
slightly negative, but a careful examination shows that 
the height of the dip begins to saturate starting around 
L = 28. In addition to the depth of the negative region 
being bounded, its width also vanishes in the thermody- 
namic limit. Finally, it is also worth pointing out that 
the behaviour of the SU(A^) transition reported in this 
work bears little resemblance to the transition in mod- 
els whose Neel-VBS transition is known to be first-order 
(e.g., Ref. [13}, which is marked by strong hysteresis ef- 
fects even at small lattice sizes. These observations sug- 
gest to us that the behaviour is only a finite-size effect. 
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FIG. 7: The spin gap As versus the inverse linear size 1/L 
close to the critical point. The dashed line indicates a lin- 
ear fit for N — 4.6 that extrapolates very close to in the 
thermodynamic limit. 

amine the transition for signs of first-order character and 
have concluded that it is almost certainly continuous. We 
have generated energy and order parameter histograms 
with extremely good statistics close to the transition and 
have found no evidence of a bimodal character. See, e.g., 
the distribution of the staggered magnetization, shown 
in the bottom panel of Fig. |3] which is single-peaked and 
evolves smoothly across the transition. 

One point of concern was that the Binder cumulant has 
a small region around the transition where it drops below 
zero. This is sometimes a signature of a first-order transi- 



VI. CONCLUSION 

In conclusion, we have introduced an algorithm to sim- 
ulate Heisenberg SU(iV) models (for the representation 
with a single row and column on one sublattice and its 
conjugate on the other) on any bipartite lattice. It is for- 
mulated in the total singlet sector and allows for efficient 
computation with arbitrary, continuous values of N. For 
the square-lattice model, we find a second-order phase 
transition between a Neel and VBC columnar state at 
Nc = 4.57(5). Constructing a Ginzburg-Landau theory 
for this phase transition is not simple from the symmetry 
point of view, as the external parameter N of the SU(A^) 
symmetry is tuned artificially to unphysical values in our 
numerics. Naively, the ingr edie nts seem similar to those 
encountered in DQC point^^^ a continuous Neel- VBC 
transition, driven by an external parameter that favors 
short VBs of the VBC over the long VBs needed for mag- 
netic ordering. Strictly speaking, the arguments of Ref. |7] 
do not apply here, since they rely heavily on Berry phase 
effects specific to S = 1/2. Hence, our results are not 
directly related to those of Refs. [TIandini The large- 
techniques of Ref. |2] do predict a continuous transition 
from Neel order to disorder, but the ground state degen- 
eracy can only be computed for integer A^. We hope that 
our numerical results for the critical exponents encour- 
age others to pursue extended analytical calculations: the 
only available estimates of exponent^^^ are for represen- 
tations with n ^ 1 and do not go beyond order 1/A^. 

Finally, we note that the algorithm presented here can 
be applied with minor modifications to the case of one- 
dimensional Heisenberg models with another generalized 
symmetry, namely SU(2)fc. This open the door to the 
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numerical study of topological quantum liquids, as found 
in Ref. 16, 

During the completion of this work, a more efficient 
algorithm for the N = 2 case was proposed by Sandvik 
and Evertz."'^^ This algorithm performs nonlocal moves 
by flipping spins around the loops in a mixed spin-VB 
representation of the projection. It is straightforward to 
generalize this algorithm to the SU(A^) case for integer 
N, as has been done in Ref. HSl One of u^^^l has recently 
shown that an algorithm for real N can be constructed 
using a loop representation of the SU(A^) model matrix 



elements similar to the one presented in Ref. Early re- 
sults obtained with this algorithm are in agreement with 
those presented in this work. 

We thank K. Harada and N. Kawashima for fruit- 
ful exchanges. Some calculations were performed us- 
ing the ALPS libraries (Ref.[2l|). We thank IDRIS and 
CALMIP for allocation of CPU time. Support from the 
Procope (Egide), the French ANR program under Grant 
No. ANR-08-JCJC-0056-01 (FA, MM, and SC) and from 
the Alexander von Humboldt foundation (KSDB) is ac- 
knowledged. 



* Electronic address: |kbeach@phys.ualberta.caj 

^ A. Auerbach and D. P. Arovas, Phys. Rev. Lett. 61, 617 
(1988); D. P. Arovas and A. Auerbach, Phys. Rev. B 38, 
316 (1988). 

^ N. Read and S. Sachdev, Phys. Rev. Lett. 62, 1694 (1989); 
Nucl. Phys. B316, 609 (1989); Phys. Rev. B 42, 4568 
(1990). 

^ N. Kawashima and Y. Tanabe, Phys. Rev. Lett. 98, 057202 
(2007). 

* K. Harada, N. Kawashima and M. Troyer, Phys. Rev. Lett. 
90, 117203 (2003). 

^ N. Kawashima and K. Harada, J. Phys. Soc. Jap. 73, 1379 
(2004). 

® G. Santoro, S. Sorella, L. Guidoni, A. Parola, and E. 

Tosatti, Phys. Rev. Lett. 83, 3065 (1999). 

T. Senthil, L. Balents, S. Sachdev, A. Vishwanath, and M. 

P. A. Fisher, Science 303, 1490 (2004); Phys. Rev. B 70, 

144407 (2004); J. Phys. Soc. Jpn. Suppl. 74, 1 (2005). 
^ K. S. D. Beach and A. W. Sandvik, Nucl. Phys. B 750, 

142 (2006). 

^ L Affleck, J. Phys.: Condens. Matter 2, 405 (1990). 
° A. W. Sandvik, Phys. Rev. Lett. 95, 207203 (2005). 
^ S. Liang, B. Dougot, and P. W. Anderson, Phys. Rev. Lett. 
61, 365 (1988). 

^ K. Vollmayr, J. D. Reger, M. Scheucher, and K. Binder, 

Z. Phys. B 91, 113 (1993). 
^ K. S. D. Beach and A. W. Sandvik, Phys. Rev. Lett. 99, 



047202 (2007). 

O. L Motrunich and A. Vishwanath, Phys. Rev. B 70, 

075104 (2004); A.W. Sandvik, Phys. Rev. Lett. 98, 227202 

(2007); R.G. Melko and R. K. Kaul, Phys. Rev. Lett. 100, 

017203 (2008); F.-J. Jiang, M. Nyfeler, S. Chandrasekha- 

ran, and U.-J. Wiese, J. Stat. Mech. P02009 (2008) 

B. I. Halperin, T. C. Lubensky, and S.-K. Ma, Phys. Rev. 

Lett. 32, 292 (1974); V. Y. Irkhin, A. A. Katanin, and M. 

L Katsnelson, Phys. Rev. B 54, 11953 (1996). 

A. Feiguin, S. Trebst, A. W. W. Ludwig, M. Troyer, A. 

Kitaev, Z. Wang, and M. H. Freedman, Phys. Rev. Lett. 

98, 160409 (2007). 

A. W. Sandvik and H. G. Evertz, arXiv:0807.0682. 

J. Lou, A. W. Sandvik, N. Kawashima, arXiv:0908.0740. 

K. S. D. Beach, (unpubhshed). 

M. Aizenman and B. Nachtergaele, Comm. Math. Phys. 
164, 17 (1994). 

F. Albuquerque, F. Alet, P. Corboz, P. Dayal, A. Feiguin, 
S. Fuchs, L. Gamper, E. Gull, S. Grtler, A. Honecker, R. 
Igarashi, M. Korner, A. Kozhevnikov, A. Lauchli, S. R. 
Manmana, M. Matsumoto, I. P. McCuUoch, F. Michel, 
R. M. Noack, G. Pawlowski, L. PoUet, T. Pruschke, U. 
SchoUwock, S. Todo, S. Trebst, and M. Troyer, J. Magn. 
Magn. Mater. 310, 1187 (2007); M. Troyer, B. Ammon 
and E. Heeb, Lect. Notes Comput. Sci., 1505, 191 (1998); 
see http : / /alps . comp-phys . org. 



